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ABSTRACT 

In this paper we reconsider a series expansion for a dark matter distribution 
function in the spherically symmetric anisotropic limit. We show here that the 
expansion may be renormalized so that the series does converge in time to an 
estimate of the steady state distribution function in the central regions. Sub- 
sequently we use this distribution function to discuss the nature of the central 
equilibrium and, by invoking stationarity of Boltzmann's H function as a mea- 
sure of (thermodynamic) relaxation, we calculate the adiabatic variation in the 
local logarithmic slope of the mass density. Similarly the pseudo (phase-space) 
density variation with radius is calculated. These are compared to empirical fit- 
ting functions. There is general agreement on the inner part of the logarithmic 
slope of the density and of the inner profile of the pseudo-density power law, but 
coincident continuity with the outer power-laws is not yet achieved. Finally some 
suggestions are made regarding the actual microphysics acting during the non- 
equilibrium approach to relaxation. In particular a cascade regime is identified. 

Subject headings: cosmology: theory — dark matter — halo formation 



1. INTRODUCTION 



This paper continues the study of the dynamical relaxation of collisionless matter using 
the analytical technique summarized recently in Henriksen (2006b: paper I hereafter). The 
reader is encouraged to refer to paper I for a more detailed introduction to the method, but 
we repeat the general idea here for convenience. 

We transform from the usual phase-space variables plus time to what might be called 
'Lie' variables in that they remain invariant under a Lie displacement (see Carter and Hen- 
riksen, 1991 for a formal treatment; or Henriksen, 1997 for pragmatic application to the 
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present problem). The direction of this displacement in phase-space-time is taken to be 
along one of the variables (T in the present paper) and the symmetry is only imposed by 
insisting that Ot = 0. Otherwise the transformation is merely a transformation to rotated 
variables in phase-space-time with no loss of generality. 

In the application to the dark matter halo problem, it is convenient to remain close 
to the homothetic or self-similar symmetry . This is a Lie symmetry wherein the relation 
between the physical quantities and the invariant quantities is found essentially by dimen- 
sional analysis. Because of the inevitable presence of Newton's constant G, the complete 
set of dimensions is reduced essentially to space and time. These arc scaled independently 
by the parameters 6 and a respectively, but all physical results depend only on the ratio 
a = a/5. 

More explicitly,Lie scaling in the dimension space of time, space and mass is represented 
by the scale vector (q;,5, /x). This can be reduced to the scale vector A = {a, 5) since the 
mass scale-factor /j, = {35 — 2a) due to the invariance of Newton's constant, G. The usual 
dimensional analysis of each dynamical quantity q is effected by assigning the appropriate 
dimensional co- vector to q. Thus for a velocity d„ = (—1, 1), while for density = (—2, 0) 
after eliminating /i in favour of 35 — 2a. 

The independent variables that express this transformation are 

R = r e-("^)/^ Y = Vr/iat)^^'"-^^ 
Z = f/{atf^/''-^\ aT = \nat, (1) 

where the usual notation for radius, radial velocity, squared specific angular momentum, 
time, and the distribution function arc used. These variables are defined to be indepen- 
dent of T. The variation of any quantity q under the scaling transformation is given by 
dxq — {dq.A)q, which integrates to g = q{R, Y, Z; T) exp (d^.A T). Thus taking q — G and 
observing that da — (—2, 3, 1), we deduce that 35 — 2 a = for the invariance of G under 
the scaling as used above. 

The strict self-similarity (now identified as a Lie symmetry) is given by holding a con- 
stant and setting 9t = 0. The advantage of these variables is therefore clear in this hmit 
since the problem is now rendered 'stationary' although it remains multi-dimensional. By 

remaining 'close' to the self-similar symmetry we mean that there is only an adiabatic vari- 
ation in a locally due to dynamical modification of the initial conditions. We recall that for 
a strict self-similarity set by the cosmological boundary conditions P{k) oc /c", where n is 
appropriate to a given range of A;, one finds (e.g. Henriksen and Widrow, 1999; Hoffman and 
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Shaham, 1985) 

a = 4i±4. (2) 

where x = 4 around a local maximum (Hoffman and Shaham, 1985) or x = 5 around an na 
peak (Henriksen, 2006b). 

Except in situations of extreme phase-space symmetry this formulation does not lead 
to exact solutions of the Poisson-Boltzmann set of equations even when Ot = . To proceed 
further analytically we make use of the freedom of either one of the parameters a or 6, here 
taken to be a, to regulate the resolution in phase-space and in time according to the relation 
for the phase-space volume 

Ar A^;^ Af = AR AY AZ e^^/"-^)"^, (3) 

and the to definition of T in terms of t. By allowing a to be large we coarse-grain the 
phase-space resolution as well as degrade the temporal resolution (an averaging). Thus this 
transformation is necessarily non-canonical from which fact it derives its usefulness. 

We thereby obtain a convenient and wcll-dcfincd coarse- graining procedure, which we 
apply by writing all quantities as a series in inverse powers of a and solving progressively for 
the various orders. In the past we have terminated the series at a given order by requiring 
all higher order terms to vanish. However the approximation found in this way diverges in 
time at each radius for any finite a. Thus our first task in this paper is to find a suitable 
renormalization procedure that renders high order terms finite asymptotically in time at 
finite a. 

In paper I the distribution function of dark matter was studied in just this manner. The 
series expansion for the distribution function / takes the explicit form to second order (see 
equation (18) of paper I) 



7rf(r,vr,f;t) - A^-'\PUCi,ClT) - u Pn{Ci,ChT) + {u' /2) P22(Ci,Cl;r), (4) 
where u is defined as 

u=l/{aR"-)^t/r'', (5) 

and the factor tt is for convenience only. We use the variable u in the renormalization 
procedure below. 

This expression uses the dimensional scaling 

7r/ = P(i?,y,Z)e-(^/"-^)'^^, (6) 
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where 

P = Po + Pi/ a + Pa/a' + • • • , (7) 

so that 

= Poor("-=^\ Pi = -Pnr("-3)/P«, 

P2 = P22r("-'V(2i?'"). (8) 

The functions Pjj are functions of the new variables ([2]) in the combinations 

Ci = rp("-^\ cl = ^P'^""'^ (9) 

These latter variables are constants on the characteristics of the coUisionless Boltzmann 
equation (CBE) to all orders. 

When approximate or 'adiabatic' self-similarity is assumed, the dependence on T indi- 
cated in the functions Pjj is weak (absent for strict self-similarity). However this dependence 
will be allowed in the renormalization procedure. The parameter a is the self-similarity in- 
dex. In the case of self-similarity determined by initial conditions it gives the ratio of the 
temporal to the spatial powers in the dimensions of a governing constant. However as argued 
in paper I and in Henriksen (2006a), this dominant constant can change dynamically. It is 
then best left as a fitting parameter. We discuss its variation further in susequent sections. 

The functions Pqo, Ph and P22 are related from the expansion of the Poisson and 
Boltzmann equations according to 



Pn = {{a - ml - 27o) + C2) + 2(a - 2)CiC2%no - Ci(3 - a)P„„, (10) 

and 

P22 = ((a - l)(Ci' - 27o) + (2)%^!! + 2(a - 2)CiC2'5c|^ii - 3Ci^ii + (3a - 2)7i%P,„. (11) 
The constants occurring in the potential, 70 and 71, are defined as in paper I by 

^° = 2(a-l)73-2a)' ^^^^ 
and ^ 

^ 3(l-a)(2-3a)' ^^^^ 
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where 

= J P,, dCidCl (14) 

The potential is given to first order by 

$ = -r2(i-")(7„ + 7,w + ), (15) 

where the general dimensional scaling between $ and the 'invariant' potential \& is 

^(r,^) = ^(P,T)e2(^/"-^) (16) 

The function \1/ has been expanded as \I' = \l'o + \E'i/a + . . . and \l/o and can be identified 
from equation (fTSll . 

The series corresponding to equation (jlj) for the density is (to second order) 

p = r-'-(I^,-IuU+^-fu' ), (17) 

where the relation between p and the invariant density B is 

p = (18) 

and 6 has also been expanded as ©o + ©i/a + . . .. 

One notes however that in the context of adiabatic self-similarity where a = a{r) the 
dependence of the density on r enters through the dependence of the lij on a as well as 
through a(r) in the exponent. 

Normally we terminate the series by insisting that P22 = (paper I) whereupon equa- 
tions ffTOl) and ffTTl) become coupled partial equations for Poo and Pn. However we proceed 
differently here. 

The difficulty with the series expansion (jl]) lies in the time dependence of u. This 
quantity increases monotonically in time and so each of the terms in the series eventually 
becomes too large to be neglected. Hence one can only expect transitory validity to such an 
expansion. 

It is important to do better than this since this series was used to discuss the power-law 
behaviour of the phase space pseudo-density (paper I). It so happens that recently McDonald 
(2006, 2007) addressed this problem in the context of the cosmological perturbative approach 
to dark matter clustering. This technique replaces the different transitory estimates of the 
unknown function by an envelope to all such transitory functions. Unlike the individual 
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transitory functions, the envelope is valid asymptotically. This technique is very well de- 
scribed in the paper by Kunihiro (1995), and we adapt this method to our problem in the 
next section. This renormalization procedure is somewhat technical and does not change 
our previous conclusions. The reader may wish to skip to section (4.1) where the renormal- 
ized distribution function is given together with a brief description of the method. The last 
section of the paper together with the discussion introduce our new physical ideas and give 
numerical results. 



2. Geometrical Renormalization 

We take the approximate functions of u for which we seek the envelope to be (mq is a 
parameter) to zeroth order 

r^^-''^F = Poo-Pn{u-Uo), (19) 

and to first order 

(i7-_^(a-3)pj 

^ = = -^11 + - «°)' (20) 

where F = irf. 

The procedure that yields the envelope function (Kunihiro, 1995) consists in differenti- 
ating the right-hand side of this equation with respect to Uo and subsequently setting Uo = u 
in order to find the envelope at u. This yields evidently to zeroth order 

duPoo + Pii = 0, (21) 

and to first order 

duPii + P22/2 = 0. (22) 

These last two equations replace the previous termination that use Pu = and P22 = 
to zeroth and first order respectively. Together with equation (fTO!) to zeroth order plus 
equation (ITT!) to first order we obtain a partial differential set for Pii(m, Ci, Cl); Poo{u, Ci, Cl) 
and finally -P22(^, Ci? Cl) if carried out to first order. One seeks to choose the arbitrary 
functions appearing in the solutions to these equations so that the expression for F becomes 
asymptotically independent of u after setting u = Uo- This latter identity is the reason we 
obtain F only to zeroth and first order even though we retain expansion terms up to first 
and second order respectively. 
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Equation fl2T]) together with equation ffTOj) yields the relevant equation for the Poo and 
Pii which we write here for explicitness as: 

((a - l)(Ci - 27o) + C2) %Poo + 2(a - 2)CiC2%^oo - Ci(3 - a)Poo + duPoo = 0. (23) 

This procedure will give the original series for the DF to zeroth order in a renormalized form 
as discussed in the next section. 

If we wish to continue to first order renormalization, we must replace the final term on 
the left in equation (l23l) by —Pu to obtain 

((a - l)(Ci' - 270) + C2') %Poo + 2(a - 2)CiC2'9c|Poo - Ci(3 - a)Poo - Ai = 0. (24) 

and then add equation (!22|) combined with equation (ITT!) in the form 

2a„Pn + ((a-l)(Ci'-27o)+C2)%Ai+2(a-2)CiC2'9c|Pii-3CiPii + (3a-2)7i%Poo = 0. (25) 

After equations (124|) and fl25|) are solved for Pqq and Pn, equation (!22|) yields P22. However 
this latter term will not appear in the DF expression after setting u = Uo- 

In practice, we normally use the zeroth order solution from equation fl23|) in equation 
fl25|) to start an iteration process that yields a non-zero Pn. The process is continued by 
using in equation fl24l) the Pn found in this way to obtain a new Poo- We shall discuss these 
renormalizations in turn in the following sections. 



3. Zeroth order renormalization 



We solve equation fl23l) by the method of characteristics, so that we must solve the set 

'^^^ {a-l){Cl-2^o)+Cl 



du 

^ = 2(a-2)CiC2', (26) 

^ = (3-a)CiPoo. 
du 



These integrate to give 



Poo = K{k^,K2)£^o''^'-''\ (27) 



where 

K^^8o{Cl)-^'^^\ (28) 
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and K2 must be found from 



dC. 



±2{2-a)u = K2. (29) 



2«:i(Cl)^^^-C|-2|7o|) 
The function K{k,i, ^2) is arbitrary and 

£o^\^^^-lo\. (30) 

The integral in equation fl29l) is not in general expressible in simple terms (although it 
is in the special case a = 0). However the entire interesting range of a is [0, 3/2] (the outer 
density logarithmic slope being —3) and in this paper, in order to describe a flattened density 
profile (relative to —2), we may consider here a < 1. Then in this domain we observe that 
the integral can be expected to be periodic in C2 since the square root will impose inner and 
outer 'turning' points. The quantity ^2 is not then an isolating integral, but gives rather the 
temporal behaviour of the (scaled)' angular momentum' on the characteristics, which are 
not in these variables the particle trajectories. For our purposes the important behaviour is 
simply the monotonic linear dependence on u. 

If at this order we wish to have an asymptotic approach to the steady state, then 
remembering equation fl27|) we may require for example that 

K(,u..) = K.MM^^^^£^). (31) 

But any form of K2 and its argument that tends to a constant asymptotically will do. If in 
addition we require a steady, spherically symmetric and velocity isotropic asymptotic DF, 
we should set Ki = constant . Hence the zeroth order, that is steady state, distribution 
function becomes 

F = r^''-^)K8i^\ (32) 

which is explicitly 



^/ = Coo|| + ^ + $o(r)|^, (33) 

as in paper I with Coo constant. This is thus our best estimate for the ultimate DF in an 
isotropic, spherically symmetric, adiabatically self-similar (i.e. a varies slowly in a crossing 
time) dark-matter core. 

We do not thus change our conclusions in paper I, but we obtain some assurance that 
the series can asymptote to the steady value. 



- 9 - 



4. First order renormalization 

If we wish to find a stable approximation to the approach to the steady state, then we 
must renormalize the first order as well as the zeroth order simultaneously. This requires 
solving equation fl25|) using the zeroth order renormalized Poo- Hence we solve the set of 
characteristic equations (see also paper I); 

2{a - 2)CiC2', (34) 

3Cl-Pll - (3a - 2)7i(%Poo)cfear, 



ds 

dCl 

ds 
dPu 
ds 



where Poo is taken from the zeroth order renormalization f l271) and ds = du/2. 

As in the previous section these equations yield the characteristic constants ki and k,2 
(although with u in equation (!29l) replaced by s). The equation for Pn on the characteristics 
becomes, by using the third member of equations (IMl) . equation (127|) . and taking (2 the 
variable, 

CI ^= (35) 



dC, 



^ -^11 + §^-^7ii^('^i, K2)Si^-'^ {p + In {K) - Krd^, In {K)g{Ci)) 



2(2 -a) 2(2 -a) 



where 



and 



(3 — 0.) , , 



s(Ci) ^ [' dx — ^ — — . (37) 

J x^-''{2kix^-'^ - X - 2\^o\) ' 

It becomes clear from equations fl20|) and fl35l) ( remembering the linear, monotonic 
dependence of K2 on u) that only the choice 

K{k^,K2) =K,{k,)/k2, (38) 

will produce an envelope behaviour (P independent of u = Uo) to first order. This also 
implies that the term proportional to dt^,^K in equation fl35l) may be dropped asymptotically 
at large u. In this case the integration for Pn proceeds as in paper I (e.g. equation (35) in 
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that paper) so we obtain formally (the constant of integration on the characteristics K{k,i, ^2) 
is also written as K = Ki{ki) / K2) 



Pn = - ( liK,{K,){p + + k,{K,)£'o'--'^ . (39) 

K2 \ duiKi J 

In the special case a = 2/3, 71 must be replaced formally by — (3/ii/2) ln£^o. In any case 
we recall that 71 oc In and hence to Pn as written immediately above, which vanishes 
asymptotically with 1/^2- Thus we may ultimately neglect this term compared to the second 
term in the bracket of equation fl39p . However it is useful to retain it as a transitory term 
(see below) in the resulting distribution function. 



Strictly asymptotically then we deduce for P] 



11 

Pu = -ki{K,)s!^ . (40) 

We turn now to the modified equation fl23l) to obtain the first order correction from 

dP 

-^ = (3-a)CiPoo + Pii, (41) 
ds 

Poo = Coo{s){(i)-^, (42) 



which, with the substitution 
becomes 



^ = iClr^Pu. (43) 
To complete the explicit solution for Poo is evidently complicated since we require 

^ = 2(a - 2)OC2^ (44) 

from equations and 

Ci' = 2ft:i(C2')i^-C2'-2|7o|, (45) 
from the integral (l28l) and the definition ( |30l) . Moreover the integral ki allows us to eliminate 
So in favour of C| from Pn. The integrals over (2 not expressible in general in terms 
of elementary functions, but they are in any case periodic because of the turning points in 
(i. Ultimately they provide too much information, since Pn vanishes asymptotically. This 
means that Coo is asymptotically constant and thus Poo retains the same zeroth order form 
as in the previous section. Consequently the distribution function renormalized to first order 
is 

F = r('^-3) (CooiK,)£f^ - uPn^ . (46) 

This is translated into the physical form that we use in our numerical discussions in the next 
sub- sect ion. 
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4.1. Summary of Renormalization Result 

The renormalization procedure followed in the previous sections due to Kunihiro (1995) 
comprises first finding an approximation to some function of interest (found generally by a 
series expansion as in our equation (jlj)). In general the approximation may approach the 
desired function only in a transitory fashion in the independent variable [u in our case). 
Subsequently the approximation is parameterized as in our equations (fT9l) and (!20|) in such 
a fashion that the value of the parameter {uo in our case) determines the region over which 
the approximation holds. The envelope to this set of functions (created as Uo varies) is then 
found by standard methods of the calculus plus the ingenious device of setting Uo = u after 
differentiating with respect to this parameter. 

The procedure is somewhat special in our case in that the terms in the approximation 
(the Pjj) are related already through equations ffTOl) and ffTTl) . and the renormalization con- 
ditions fl?T]l . fl2^ merely close the system at the respective orders. The previous sections 
have thus been occupied with the solution at zeroth and first orders for P^o and Pu which 
render the corresponding series (fT9|) and (l20l) independent of u. 

This procedure concluded with the result of equation (146|) (together with equation (l39l) 
which in terms of physical quantities becomes explicitly 

ttA/ = (47) 
CooA\E\^^) - ki{Ki)\E\^(^)r'' - -tiKi{Ki) (^^Jriy + \E\^^) 



l^l = ly + ^ + <fo(r)|, (48) 
a ± sign has been absorbed in the constants and A = 2(2 — a). 

Once again for a = 2/3, we must formally replace 71 by — (3/ii/2)(ln + lnr~^/^), 
although this term vanishes in time in any case as remarked previously and is strictly of 
transitory interest at a fixed r. 

Equation (H7|) gives the formal result of this paper. Thus asymptotically, when the time 
dependent term disappears, the first order term does not dominate the zeroth order term as 
r — > 0. This was not the case for the original unrenormalized series (jll). There is however 
still a limit to the validity of the asymptotic (in time) expression at large r where the first 
order term becomes comparable to the zeroth order term. 

We have left the decreasing time dependent term in 71 (actually In) because of the r 
dependence at a < 1. This may be used to off-set the time dependence (that is we need 
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j,(2-2a) ^ which is r oc t^/^^"") to keep this term's importance relative to the zeroth order). 
Thus this term allows some estimate of the approach to the asymptotic form at any fixed 
r. The dominant first order term proportional to Ki is renormalized in time but as already 
noted gives an outer limit to the validity of the series when it becomes comparable to the 
zeroth order. We normally take this term to be negligible at small r in what follows and 
give our preliminary applications uniquely in terms of the zeroth order (equilibrium) part of 
( 147|) . The time it takes for this equilibrium to be attained may be estimated in principle at 
each r by insisting that the transitory term be negligible, but this would require a precise 
determination of the constants for a particular system. 

In the next section we discuss especially the implications of the zeroth order for the 
central density and pseudo- density profiles, including the adiabatic variation of a. 



5. Adiabatic Approach to Central Equilbrium 

Taking the zeroth order result from the previous section as the central distribution 
function of a dark matter halo, we calculate the density from the integral 

p(a,r) = 2^/l2)Coo / dE ^{E - $„) (49) 

which yields after the substitution u = ^o/E 

p{a,r) = 2v^Coo|7o|^r-2'^ D{a;Um), (50) 
where the ad hoc function D is 

/ . 2a-l 

D{a;Um) = / duy/l — uu^-^, (51) 

Dia; u^) = Bi-^, 3/2) ut~^ 2^^1(7^, -1/2; ; u^). (52) 

I — a a 1 — a 1 — a 

We define 

Um = ^o/E^, (53) 
where Em is the maximum energy present at r. 

The quantity |7o| is found by completing the integral for I 00 as in paper I and using the 
definition (fT2|) to be 

'^°'^-° = (l-a)(3-2a) ^^^^ 
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and this in turn gives ^o{o,,r;Um) from equation f|T5l) . 

In these expressions B{x, y) is the Euler Beta function and 2-^1 is the Gaussian 

hypergeometric function. 

This gives the density exphcitly as a function of a and r in the form 

p = {kDY^''\{l - a)(3 - 2a)y r'^^ (55) 

Hence the logarithmic slope j3 = —dlogp/dlogr is found as 

(3 =2ax 

1+— In A;D a,f/™ r In 3-5a+3a )+ 2^ 

2 eta (3 — 5a+2a^) 

for a < 1, where we have defined ax = dlna/dlnr and k = 2\plCoo- 

Clearly if is small /3 ~ 2a(r). This tends to be the case except at small /c, which 
correponds to a small central density and/or at relatively large r. In such cases our argument 
below yielding air) in terms of the Boltzmann H function, probably fails in any case because 
of insufficient relaxation. 

Indeed an important consideration is that only with Um 7^ does one obtain a finite 
value for p as a,r both tend to zero. Otherwise the Beta function in -D(a, Um) dominates and 
the density diverges as 1/a {B{x, 3/2) ^ 1/x as x — > 0). Although with Umc = lirrir^oUm = 
const. 7^ the central density is always finite, it becomes only logarithmically larger as m^c 
approaches zero since the cancellation of the 1/a divergence in D requires |alnMmc| ~ 1- 

However, using <l>o in the definition explicitly, Umc = hm(^^o)(|7ok^''^~'*V-^m) where Em 
is the maximum energy present at r. Unless Em — in the same limit Umc goes to zero as 
r — >■ (a < 1), which by the above discussion would imply an infinite central density. We 
therefore conclude that Em oc $0 as r — >■ (so that both quantities vanish at the centre at 
the same rate) if one is to have a finite central density. The energy dispersion (and hence 
the velocity dispersion) of the central particles will thus decrease ultimately towards zero. 

From the empirical fitting formula of Navarro et al. (2004) we deduce that 

P = 2{—r-, (57) 

where r_2 is the radius where /? = — 2 and lip 0.17, although there is dispersion about 
this value at the 20% level. If one ignores the terms multiplying a-p in equation fl56p we find 

« = (— r- (58) 
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Thus a and r do tend to zero together (and a is adiabatically slow relative to r) according 
to the numerical simulations and consequently in our expression Um 7^ 0, however small, 
in order to guarantee a finite central density. One can imagine that such delicate physical 
behaviour may well be refiected in the difficult convergence of numerical simulations. 



Also from (H7!) we may calculate the rms isotropic velocity dispersion cr^ from f ^ — f ^ 
that is 



<y' = - fmE - $ o(r)))3/2 _ 2(£; _ $^)) dE, (59) 



which yields at equilibrium (zeroth order) 

2 ^ / k \ 2{i-a) I Gbia;u^) _ 2Gi{a,Un,y 



V(l -a)(3-2a)y \D{a;Um)'' Dia;UmY^+-^ ' ' ^ ' 

Here we have defined new functions 

G,{a-u^) ^ (i)[x^ ^F,{q,-^/2--^,x)t^, (61) 
q 1 — a 



where 5 = 1_ and 



^ — , (62) 



plus 

(1 + a) \[Sa - 1) J 

With the adiabatic variation a(r) the r dependence of this expression for o"^ is no longer 
obvious. Despite appearances, it is well behaved at a = 1/2. 

The quantity Umc determines the rate at which the particle energy decreases to zero at 
the centre of the system since $0 vanishes there as discussed previously. For simplicity we 
have taken Um to be constant everywhere in our numerical results rather than just at the 
centre of the system, but it might vary with r and a(r) through the dependence 

^ ^JE^ = ( (64) 

' \il - a){3 - 2a) J Em ^ ' 

if Em is instead taken constant for all r. 

As a check on this assumption we plot in figure ([T]) Mm(fl, ^) over the domain of interest 
for our reference parameters. We see that it is approximately constant at the value 10^^^ 
that we use below. The right panel shows the surface of the Boltzmann H function with 
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variable Um to be compared to the corresponding plot in figure ([2]) for constant Um- We see 
that, but for the the unphysical region at large a and small r, the correspondence is good. 
We use the Boltzmann function in our analytic arguments of the next section. 

As we discuss somewhat further in the next section, both of the parameters k, viffi are 
important to the behaviour of the density and pseudo-density. This implies a dependence 
on the central density of the system {k) and on the rate at which high energy particles are 
absent from central regions {um)- There is then strictly a lack of universality in the predicted 
profiles. However the dependence on the central density could be scaled away at least at 
small a, and the dependence on Um is relatively weak except at values much closer to unity. 

We have thus chosen a particular physical system as our reference, which has parameter 
values k = 250 and Um = 10"^*^. Such a system (see results below) yield plausible radial 
profiles (by which we mean similar to those found in simulations) both for p{r) and the 
pseudo-density 0(r) out to the radial limits of our assumptions and close to the inner limits 
of the simulations. Our outer limits are set by the efficiency of the relaxation and to extend 
them will require adding in the higher order terms in equation (14 7p . 




Fig. 1. — The left panel shows the surface Um{(i, r) over our domain of interest with k = 250 
and Em = 10^^. The right panel shows for the same parameters the surface of constant 
Boltzmann H function to be compared to the corresponding plot in figure ([2]). 
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The 'pseudo phase-space density' (j) (Taylor and Navarro, 2001; Hansen, 2004) may be 
found simply from equations fl55l) and fl60l) as (strictly e = 1 for their quantity) 

= p/a^\ (65) 

The explicit dependence on radius is oc (e.g. paper I and Hansen (2004)) 

n = (3e-2)a-3e. (66) 

As for the density the variation of (j) with r is no longer simply given by the power n because 
of the locus a(r). 

However from the analytical treatment one would like to deduce the variation of (3 
or even a, at least over the limited range in r wherein the fitting formula is verified. To 
this end we examine the Boltzmann H function (an entropy density) that follows from our 
zeroth order equilibrium. We might expect the microphysics of dark-matter halos to produce 
relative constancy of this function along the locus a(r), since a is monotonic in time at each 
r until equilibrium is reachecl^. 



5.1. H Function and the Locus a{r) 

Although it is understood that a self-gravitating system may not have a condition of 
maximum entropy globally since this leads to an isothermal DF and an infinite mass, it 
is possible that locally (e.g. in the central regions) a dark matter halo may maximize an 
entropy density. This leads us to consider the stationary behaviour of the Boltzmann H 
function that follows from the zeroth order of (H7|) . This is formally 

H, = j foinfodh. (67) 

Proceeding with a straightforward but tedious calculation of the integrals we obtain explicitly 

Ho{a,r;Um)= (68) 
AV2Coo $1"^^ ((In ^ + 1^ In <^o)D{a; + £^^^^9, D{a; uj) . 

This surface is shown in figure ([2]) over an interesting range in a, and r. There is in 
addition an arbitrary fiducial unit for r. The parameters Um = 10~^^ and k = 2\piCoo = 250 



recent paper (salvador-Sole at al., 2007) suggests that essentially equation ([2]) with the appropriate 
history of n is sufficent to determine a(r). This would oppose the notion of universal central relaxation used 
here, although the two assumptions appear to yield compatible results in the core. 
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as noted above. We begin to see the divergence at small a that exists strictly for zero m. 
The accompanying graph in figure ([2]) shows the behaviour of the function F for the same 
Um (independent of k). The limiting value as a ^ is ~ 40.833 

We observe that, by using the equilibrium distribution function that is the zeroth order 
of result (H71) . we calculate an Ho that is approximately constant compared to the value at 
large a and small r. This latter case does not arise physically. The behaviour is what one 
might expect in an equilibrium state. Our additional hypothesis is that we may deduce the 
variation a(r) by seeking the precise stationary locus SHo{a,r) = 0, for which we must solve 

^ = -M£. (69) 
da drlio 

The solution to this differential equation is shown in figure ([2]) for Um = 10~^^ in the 
left bottom panel. In this example we start at r = .01 and a = 0.5. Some variation of the 
starting values is possible. We observe that there is a natural outer limit to the usefulness 
of our procedure where a becomes double-valued. We regard this point as the outer limit to 
the necessary relaxation, since lower densities (proportional to k) causes this limit to move 
inward. This particular curve is illustrative only, and some fine tuning in k (proportional 
to the central density) and starting value would be necessary for any particular halo. The 
physical scale depends on the fiducial unit, but the relative scale is fixed. However this 
theoretical curve does indicate that a varies slowly with r, as is suggested by the numerical 
simulations (e.g. Navarro et al., 2004). 

The question now arises as to whether there is an approximate analytic representation 
of this theoretical locus. Since at small Um there is a considerable range of a over which the 
Beta function term in D{a; Um) dominates; we can hope to find an analytic expression for the 
locus in the small a region by letting a — >^ and holding Ho constant, while retaining only 
this term in D. Proceeding in this fashion, but keeping only the largest terms ( including 
large constants; as well as large constants multiplied by a, which is after all finite in the 
range of interest), we obtain an estimate of the locus r(a) as 

r = v^exp (-l/(2a))C(a), (70) 

where the variable constant C{a) is actually 

The first two factors in C are weakly dependent on k so that the constant depends on 
aHo/ {Qk) = h a.s a —>■ but finite. We show this variation in the left upper panel of 
figure ([3]) for k = 250 and we observe that but for the region of large a and small r which 
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is unimportant physically, it is relatively constant. The right panel shows the resultant 
illustrative fit {h = —1.0) to the numerical locus (solid line) as a dash-dotted line over the 
region where we believe our arguments to be justified. A locus of the type suggested by the 
simulations (Navarro et al. 2004) (although outside their stated parameter range by about 
two sigma) is shown as the dashed line. The approximation diverges from the exact locus 
of the stationary H function (solid line) where expected, but interestingly it continues to fit 
the simulation- type curve over a wider range in r. 

If one solves equation (!70l) allowing the exponential factor to vary, one finds two branches. 
One of which (the smaller r at a given a) is the branch that matches the other estimates, 
and the other is clearly unrealistic. 

In the bottom row of figure ([3]) we show on the left panel the density from equation flS^ . 
The outer limit to the assumption of equilibrium is clear where it becomes double valued. 
The slope is shown on the right hand lower panel to a ~ .35. At a = 0.4 and at a = 0.5 
the slope is 0.794 and 0.98 respectively. Thus we see that in this interesting range 2a is in 
fact a reasonable estimate of /3. Moreover the slope begins near the value 1, which occurs 
frequently in the simulations. However it rolls over to a fiatter value rather faster at first 
than seems to be indicated by the simulations. 

We should remember in considering these figures that there are arbitrary units, which 
must presumably be chosen so as to fit smoothly to the outer well-established NFW be- 
haviour. The dynamical variation of a in this region remains mysterious, although it must 
have something to do with the evolution of the radial orbit instability (see discussion). 

This rough continuity with the simulations does not hold universally however. For lower 
central density the results are dramatic. For /c = 100 the logarithmic slope never gets larger 
than 0.38 at a ~ 0.3, and 2a does not approximate (3 well until a < 0.2. Hence for 
lower densities the domain of applicability of this argument from relaxation moves distinctly 
to smaller radii and smaller j3. We might therefore conclude that the central density must 
be sufficiently high for the relaxation to overlap substantially with decreasing density. This 
seems intuitively reasonable, but does not seem to be sufficient to guarantee the simultaneous 
presence of the power-law behaviour in the pseudo-density,^ . 

This puzzle arises as we consider both the isotropic velocity dispersion cr^ and 'pseudo 
phase-space density' 0. Plots of these quantities are shown in our final figure (jl]) of this 
section, together with the linear relation between J and a and the previous calculation of a 
versus log{r). This latter pair is simply a convenient device for relating J to a, (3, and the 
spatial scale r. 

We see that the desirable power law in the pseudo-density exists only out to r ^ 10^'^ 
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where the velocity dispersion peaks. Inside this radius the density is rather flat so that 
it seems that a flat core must be forming before the pseudo-density becomes a power law 
contrary to what is observed in the simulations. 

This behaviour is characteristic over a broad range of parameter space {k and Um)- For 
example at k = 150 the pseudo-density power-law has moved out only to ~ 10~^'^ where the 
density is already quite flat. By taking k = 10 and Um = 10"'' one can move the pseudo- 
density powerlaw out towards continuity with the NFW regime, but now the density is flat 
and rising slightly (not really a power law but a typical logarithmic slope is —0.2). This may 
be an effect of holding Um constant at too large an r, since effectively this cuts off the high 
energy particles( We do expect this to happen as r — > to ensure a finite density). Indeed 
one finds that holding constant leads in this case to a substantial variation (factor 4) in 
Um in our domain of interest at finite r. 

Clearly both k and Um are significant parameters and more detailed study is required. 
Such a study would make the present work far too long. It probably will not be solved by 
including the first order term in fH7|l . since this was done in paper I with similar results. For 
the moment it seems that either the density in the simulations should flatten rapidly in the 
next order of magnitude in decreasing r, or the power-law behaviour of the pseudo-density 
should break in the same region; if our hypothesis is correct and/or the measures from the 
simulations are consistent. 

We conclude from this preliminary study that at least qualitatively the stationarity of 
the H function seems in accord with the empirical behaviour found in the simulations at the 
inner limit. However the details are discordant in general and there is a lack of universality at 
different central densities, although this is mainly a question of changing the scale at which 
similar behaviour appears. The microphysics of the inferred relaxation remains somewhat 
mysterious, although we believe it to be intimately related to a cascade from the radial orbit 
instability. We discuss these points further in the concluding section 
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6. Discussion and Conclusions 

A key technical result of this paper is a successful renormalization in u of the series 
that underlies the basic calculation of paper I in section 4.2. This reassures us that the 
consequences of using a divergent series to study the pseudo-density has not led us astray. 
For the critical quantity termed 'ratio' in paper I is, asymptotically in time, independent of 
the renormalization factor (i.e. 1/^2). The conclusions drawn there are thus qualitatively 
unchanged by the renormalization. Even though in that paper we had no natural way of 
calculating a(r) the results at a given a are rather similar to those found here. 

We have also confirmed the asymptotically time- independent form (zeroth order DF) 
for the distribution function at fixed r. This is defined to within arbitrary functions of ki 
in principle, but the requirement of velocity space isotropy renders these as true constants. 
The development of isotropy requires a physical mechanism for expanding the phase space 
available to the system. This we attribute to the radial orbit instability (see e.g. MacMillan, 
Widrow and Henriksen, 2006). This equilibrium distribution function is stable according to 
the usual theorems (e.g. Binney and Tremaine, 1987). 

There is a first order correction to the strict equilibrium DF which is now renormalized 
in time. This permits an estimate of the size of the relaxed core (a 7^ 0) by taking it equal 
to the zeroth order. The estimate is energy dependent being roughly r^ore 
although we should remember the adiabatic dependence of a on r. This result indicates that 
at any r it is the high energy particles (recall that E > for the equilibrium result) that are 
the more likely to be relaxed. This is consistent with such particles having large turn-round 
radii even if they are relaxed near the centre of the system. 

Our second result is the deduction from equation (!69l) of a form for the variation of the 
logarithmic slope of the density (related to a(r)) that is not qualitatively different from that 
found empirically (Navarro et al., 2004). The expected qualitative behaviour with r is also 
inferred for the velocity dispersion and the pseudo phase-space density. 

However there is a slight contradiction between this latter behaviour and that of the 
density. The pseudo-density is more like the simulations at a given central density for lower 
energy particles at the centre (large Um), while the density profile is not so dependent on 
this parameter. There is some evidence that there may be an optimum value for matching 
the density profile to the simulations however [um = 10~^^ is slightly better than larger and 
smaller values), which is in accord with the lower energy particles being more relaxed. 

We have concluded from these considerations that, assuming correctness on both the 
theoretical and simulative sides, either the density profile or the pseudo-density profile should 
change dramatically on proceeding to a scale smaller by an order of magnitude than accessible 
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currently. 



If the disagreement with the simulations persists as strongly as at present, we might 
suspect that something is missing in the theoretical formulation. It is possible that some 
form of collision term is necessary in the Boltzmann equation in order to describe properly 
the wave-particle interaction. 

We have used the stationarity of the Boltzmann H function to deduce these profiles. 
However this argument does not by itself provide a physical mechanism for the onset of such 
a condition. We must identify some scattering mechanism that increases the phase-space 
volume accessible to the system, relative to that available to nearly radial infall. 

The likely candidate appears to be the radial orbit instability (e.g. Barnes et al. 2005, 
2006; Lu et al.,2006; MacMillan, Widrow and Henriksen, 2006) plus a subsequent cascade 
to smaller scales. Such a cascade might be due to the wave-particle interaction as discussed 
in Henriksen (2006a) and elaborated below. The self-similar form of the squared angular 
momentum that follows from equation ([2]) namely = (^|r*^^~^"'' already requires something 
like the radial orbit instability to generate this quantity. Moreover we see from equation 
([71]) below that j is proportional to V Gmr, which behaviour has been associated with the 
radial orbit instability and indeed previously with the existence of self-similarity (MacMillan, 
Widrow and Henriksen, 2006 and references therein). 

The process is presumably working during the evolution from the density profile to 
that of and flatter. However our argument from 'thermodynamic' equilibrium does not 
extend to this region. The evolution must then be occurring without the local maximization 
of entropy such as represented by a cascade. 

Since the locus a(r) presumably contains the relaxation physics, it is of interest to see 
how it appears in various dominant physical quantities during the approach to the 'relaxed' 
region from the outside. 

Let us first recall that the mass m(r, t) of the growing halo is given by (see e.g. Henriksen 
and Widrow, 1999 ) 




(72) 



which on recalling equation ([2]) becomes 



m(r, t) 



M{R){r/R)^^-^''\ 



(73) 



Consequently we obtain 



Gm 



GM 



(l-2a) 



(74) 



■r 
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The first factor on the right-hand side of this last result is constant at a fixed fraction 
of the growing halo (i.e. constant R), so that the physical gravitational acceleration acting 
on particles there varies as r^~^". 

Starting at a = 1 which applies outside the isotropic core where the velocity is pre- 
dominantly radial and lies outside our calculation from the Boltzmann H function, we can 
describe a systematic decrease in the dependence of the radial acceleration on r as a de- 
creases. This should correspond to decreasing susceptibility to the radial orbit instability 
(see e.g. MacMillan, Widrow and Henrikscn, 2006). In fact however it is probably the action 
of the radial orbit instability that decreases a initially. Wc recall that by the approximate 
self-similarity we expect /3 ~ 2a in this regime with a locally constant. 

At a = 1 the radial gravitational acceleration at a constant fraction of the growing 
core varies as = R^^{at)^^. Hence it is very strong initially {t 0) at all R and 
creates predominantly radial infall with a density profile p oc r~^, a constant a, and hence a 
pseudodensity power n — —2. This state is highly susceptible to the radial orbit instability 
(MacMillan, Widrow and Henriksen, 2006). We ignore the r"^ outer cut-off here, which has 
been described elsewhere as due either to mass exhaustion or tidal truncation (r~^ in that 
case: e.g. Henriksen, 2004). 

We suppose a to decrease with r as a consequence of the instability until it reaches 
3/4, which is a hmiting value found by Moore et al. (1999). This yields a weaker radial 
acceleration proportional to r~^/^ = R~^/'^{at)~'^/^. At any fixed core fraction we should 
expect this phase at a sufficiently late time, but it arises at small core fractions first. The 
weaker radial acceleration presumably coincides with the trend towards isotropy in the core, 
since the radial orbit instability has already acted most efficiently. 

For a — 3/4 one finds also p oc r~^/^, a oc r^/^ and n — — (3/2)(l -|- e/2), by treating 
a as locally constant in the full expressions. This latter value is —2.25 for e = 1. This 
may seem too large compared to the value of n found in cosmological simulations (e.g. 
Dehnen and McLaughlin, 2005), but interestingly it corresponds quite closely to the value 
found by MacMillan (2006) in a study of the collapse of an isolated halo with cosmological 
initial perturbations. He finds globally n = —2.19 ± .03, and a logarithmic density slope 
close to P — —1.5 is a reasonable mean value over about a factor 4 in radius. Moreover a 
increases roughly as r^/^ in these core regions according to MacMillan's simulations. It is 
possible that the steeper value of n is an effect of the isolated nature of this halo compared 
to the cosmological simulation, perhaps due to a resulting lower velocity dispersion. In our 
calculations we find that (f) is steeper when is larger, that is when fewer high energy 
particles are retained. 
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Should a attain the value 2/3 by continued relaxation with decreasing r, we obtain 
locally the Evans and Collett (1997) solution, wherein the radial acceleration at a fixed halo 
fraction is cx r~^/^i?~^/^(at)^^/^. This phase again becomes more important at smaller R 
and late times. We expect the trend towards isotropy to be even more strongly developed 
here. One finds consistently for this phase with locally constant a; p oc r^'^/'^, a oc r^/^ and 
n = 4/3 + e, that is n = 7/3 for e = 1. This re-expresses the transitory steepening in n as 
described in paper I, but for which there is little evidence in the cosmological simulations. 

The intriguing part of this steady solution is that it can be interpreted as solving the full 
collisional Boltzmann equation (Evans and Collett, 1997). It is unlikely that particle-particle 
collisions play a role in the evolution of dark matter halos, but wave-particle collective 
interactions, equivalently clump- clump and sub-clump-clump interactions, may well do. 
Indeed there is a Kolmogorov cascade from larger to smaller radii (cr^/r = const) in this 
solution, which implies inter-scale coupling, and fits well with the idea that a cascade from the 
radial orbit instability is driving the evolution exterior to the relaxed core. Indeed Henriksen 
(2006a) has suggested that such a cascade is essential to the relaxation of collisionlcss matter, 
where the 'clumps' on a given scale were essentially identified with Landau-damped waves. 
Such a mechanism allied with the radial orbit instability at large scales would create the 
approach to the stationarity of the H function on this view. 

Pursuing this cascade idea slightly further, we note that a^/r oc r*^^^^'*^ locally. Thus for 
a < 2/3 the energy flux decreases to smaller r and increases to larger r while the converse is 
true for a > 2/3. 

This cascade behaviour suggests to us that a state with a < 2/3 will trap energy, 
perhaps creating an inverse cascade back to larger radii. This may increase the effective a 
(thus reducing the radial dependence in a) at larger r. But if it increases above a = 2/3, the 
cascade is now capped at larger radii and it should again reverse. In this way a = 2/3 may be 
a stable point in the adiabatic self-similar history of the core, about which slope oscillations 
may occur. This state is ultimately transitory, since if entropy is steadily increasing locally 
through the cascade process, we expect that considerations such as those given above will 
apply. This leads to a stationary H function and the decline in a will continue. 

In order to obtain the NFW profile we must consider that relaxation continues until a — 
1/2 and beyond . This limit is distinguished by having a constant gravitational acceleration 
at a fixed fraction of the growing core. It is just at the boundary of what we deem to be the 
relaxed region so that subsequent behaviour would be as discussed in the previous section. 

Overall we have thus explored several new ideas in this work. The technical part uses 
the renormalization idea pioneered by Kunihiro and McDonald (ibid). We have in addition 
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examined for the first time tfie consequences of maximizing entropy locally, given a form for 
the distribution function that can vary adiabatically. Moreover we have discussed the outer 
region where this is not likely to be true in terms of an energy cascade from macroscopic 
instabilities and the strength of the radial acceleration. The instabilities may be internal 
such as the radial orbit instability, or they may be created by mergers. The results are 
promising for our understanding, but not yet definitive. 

This work was supported in part by the Canadian Natural Science and Engineering 
Research Council. The author thanks an anonymous referee for a careful reading of the 
manuscript and constructive criticism. 
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Fig. 2. — On the left we show part of the surface of the Boltzmann H function over a and 
r for k = 250 and Um = 10~^^. The locus of interest is given by 6H{a,r) = 0. The right- 
hand panel shows F{a,Um) = D{a,Um) going to a finite limit as a ^ 0. The left bottom 
panel shows the locus a(r) of stationary values of Ho that follows from equation (169|) in 
an interesting range of a. We used k = 250, Um = 10"^*^, and r = .01 at a = 0.5 in the 
calculation. 
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Fig. 3. — The left upper panel shows the suface h = aH/k for k = 250 and Um = 10^^^. 
Taking the approximation of equation (1701) with h = —1.0 one obtains the dash-dotted line on 
the right-hand panel, while the numerical solution (started again at a = 0.5 with r = 0.01) 
is shown as the solid line. The dashed line on this panel is a = O.Tr*^'^^. The lower left 
panel shows the actual density for this case calculated from equation ([55]) as a function of 
log r. The lower right panel estimates the logarithmic slope /3 as a function of J, where 
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Fig. 4. — The left upper panel shows the variation in o"^ with r on a log-log scale while the 
upper right panel shows the 'pseudo- density' on a similar scale. The lower left panel plots 
the linear relation between the naively predicted logarithmic slope 2a and the index J. The 
lower right panel repeats the solution to equation ( l69l) for convenience. All of these graphs 
are for our reference example with k = 250 and Um = 10~^®. 



